Forecasting the combined effects of future climate and land use change on the suitable habitat of Davidia involucrata Baill

Abstract Accurately predicting the future distribution of species is crucial for understanding how species will response to global environmental change and for evaluating the effectiveness of current protected areas (PAs). Here, we assessed the effect of climate and land use change on the projected suitable habitats of Davidia involucrata Baill under different future scenarios using the following two types of models: (a) only climate covariates (climate SDMs) and (b) climate and land use covariates (full SDMs). We found that full SDMs perform significantly better than climate SDMs in terms of both AUC (p < .001) and TSS (p < .001) and also projected more suitable habitat than climate SDMs both in the whole study area and in its current suitable range, although D. involucrate is predicted to loss at least 26.96% of its suitable area under all future scenarios. Similarly, we found that these range contractions projected by climate SDMs would negate the effectiveness of current PAs to a greater extent relative to full SDMs. These results suggest that although D. involucrate is extremely vulnerability to future climate change, conservation intervention to manage habitat may be an effective option to offset some of the negative effects of a changing climate on D. involucrate and can improve the effectiveness of current PAs. Overall, this study highlights the necessity of integrating climate and land use change to project the future distribution of species.

change and has important implications for guiding future conservation planning.
Currently, it has been evidenced that global climate change and land use change are two of most important contemporary factors in driving species redistribution (Newbold, 2018;Powers & Jetz, 2019). However, previous studies focused on modeling spatially explicit patterns of species' range shifts under future environmental conditions by using species distribution models (SDMs) (Dawson et al., 2011;Pearson & Dawson, 2003), a powerful tool for projecting future species ranges (Elith & Leathwick, 2009), have so far mostly considered single factors, focusing extensively on climate change and neglecting other important pressures, such as land use change . This approach of single factor to estimate future species ranges may lead to unreliable projections (Marshall et al., 2018;Sirami et al., 2017;Titeux et al., 2017) and has recently been increasing questioned (Marshall et al., 2018;Regos et al., 2016).
In this context, several studies have attempted to project species future distribution by integrating climate change and other drivers, such as land-use change (Elsen et al., 2020;Pant et al., 2021;Zhang et al., 2017) and habitat connectivity (Brun et al., 2020;Huang et al., 2020;Yesuf et al., 2021). This type of integrated approach largely contributes to a deeper understanding of how interactions among drivers may affect the future distribution of species and accordingly is essential to set effective conservation policy and management (Di Febbraro et al., 2019;Zamora-Gutierrez et al., 2018), and thus has been widely recommended to modeling species future distribution .
As a tertiary relict plant endemic to China (Fu & Jin, 1992), the Davidia involucrata Baill. currently ranges approximately from 98 to 110°E, 26 to 32°N in southwestern and south-central China (Li, 1954;Liu et al., 2019;Takhtajan, 1980;Tang et al., 2017). It offers an especially rich system for exploring the interactions of climate change and land use change for several reasons: (a) as a rare species listed in the China Plant Red Data Book under first-grade state protection, as well as an "Endangered" species on the IUCN Red List of Threatened Species, it is important ecologically and economically (Liu et al., 2019); (b) long-term national observation records of D. involucrate are available in China (Tang et al., 2017;Wang et al., 2019); (c) previous studies have explored the vulnerability of D. involucrate to future climate change (Long et al., 2021b;Tang et al., 2017;Wang et al., 2019); (d) despite its populations are often found in subtropical evergreen broad-leaved forests or in mixed forests of temperate deciduous broad-leaved trees at high altitudes of between 1100 and 2600 m (He et al., 2004), land use change has led to a sharp decrease of its remaining habitats (Wang et al., 2019); and (e) future climate and land use data under different scenarios are available in China (Fick & Hijmans, 2017;Li et al., 2016).
Here, we aim to (a) examine the combined impacts of climate and land use change on the future distributions of D. involucrate and (b) evaluate the effectiveness of current PA networks in protecting D.
involucrate under the combination of climate change and land use change scenarios. To do so, we compared projected distributions produced by two types of models: (i) including only climate variables (climate SDMs) and (ii) adding also land use variables (full SDMs). We expect that the inclusion of land use covariates will significantly improve the performance of SDMs and that the differences in the projected distributions produced by climate SDMs and full SDMs will be significantly. By doing so, we provide insights into the conservation of D. involucrate in the face of pervasive environmental change.
2 | MATERIAL S AND ME THODS 2.1 | Study area and species occurrence data As D. involucrate is a Tertiary relict plant endemic to China, we chose the whole China as the study area following previous studies ( Figure 1) (Long et al., 2021b;Tang et al., 2017;Wang et al., 2019).
Totally, 337 occurrence records of D. involucrate at the spatial resolution of 10 km × 10 km for the time period 1979-2013 were obtained from Long et al. (2021a), in which an extensive database of occurrence records was assembled. Additional details on the methods of collecting and processing the occurrence data can be found in Long et al. (2021b). Following Long et al. (2021b), we used six bioclimatic variables at a 10 km resolution averaged for the period 1979-2013 from the CHELSA database (Karger et al., 2017) for model training and projection of the current potential distribution of D. involucrate: annual mean temperature (BIO1), isothermality (BIO3), temperature annual range (BIO7), precipitation of the driest month (BIO14), precipitation seasonality (BIO15), and precipitation of the warmest quarter (BIO18), for these variables have low multicollinearity (all variables with Pearson correlation coefficients |r| < .7 and variance inflation factor VIF < 5) and have the greatest ecological relevance to D. involucrate (Long et al., 2021b). Similarly, the same six bioclimatic variables at a 2.5 arc-min resolution for two future time periods, 2050s (averaged for 2041-2060) and 2070s (averaged for 2061-2080), under two representative concentration pathways (RCPs) scenarios, RCP2.6 and RCP8.5, from six global circulation models (GCMs):
The current and future land use data at a 1 km resolution were obtained from the FROM-GLC database (Li et al., 2016), which includes the proportion of ten different land-use types: (a) bareland,  For these algorithms require species absence records, we generated 10,000 pseudo-absence points by randomly sampling without replacement. We then performed cross-validation on each algorithm, where random subsets of 70% of data were used for model calibration and the remaining 30% for evaluating model performance. At each cross-validation step, two metrics were used to assess the predictive performance of each algorithm: the area under the relative operating characteristic curve (AUC) and the true skill statistic (TSS).

| Species distribution models
To avoid random bias, this cross-validation process was repeated 10 times.
To produce robust forecasts of the distribution of D. involucrate, we produced an ensemble model based on the algorithms with AUC > 0.80 and TSS > 0.60, which creating a median representation of the predictions of the 10 runs and the selected algorithms together. For each covariate included in the ensemble models, we also estimated variable contribution by calculating the change in correlation between the covariates and the response with and without the selected variable (Thuiller et al., 2021). The final ensemble models were then projected to the environmental conditions under current and future period, respectively. One probability map for habitat suitability was produced for each of the two model types at 2050s and 2070s under 12 future global change scenarios, representing all combinations of the two RCP scenarios and the six GCMs. Finally, all these maps were converted to binary presence absence maps by choosing the probability threshold that maximized the TSS value (Thuiller et al., 2021).

| Statistical analysis
Analyses were conducted on the above binary presence/absence habitat suitability maps. To explore the effect of climate and land use change on the projected distribution of D. involucrate, following Zhang et al. (2015), we first calculated two metrics of species' vulnerability to future global change: the relative changes in suitable area of D. involucrate between current and future periods in the whole study area and in its current suitable range, respectively, following the Equation C = (A − B)/B, where C is the relative change F I G U R E 1 Potential suitable (yellow) and unsuitable (gray) habitat suitability of Davidia involucrata Baill. In China projected by climate SDMs (a) and full SDMs (b). Red points represent occurrence records of D. involucrata ratio, while A and B are the future and current potential suitable habitat area (km 2 ) of D. involucrate. The first metric measures to what degree the future suitable area is large or smaller than the current suitable area, while the second metric measures the loss of current suitable habitat (Thuiller et al., 2011(Thuiller et al., , 2014. Second, to examine the distance and direction of spatial shifts, we also took the centroids of the species range from the current and the future periods using the R package "rgeos" with the "gCentroid" function. A positive distance value indicates northerly shift and negative, a southerly shift. Finally, to explore the conservation effectiveness of current nature reserve networks in protecting D. involucrate under global change, we obtained spatial polygon data for nature reserve boundaries in China from the World Database on Protected Areas (UNEP-WCMC, 2017), which was resampled to 1 km spatial resolution. We further calculated the relative changes in potential suitable habitat area inside the current nature reserve networks between current and future periods both in the whole study area and in its current suitable range, respectively.

| Species distribution models
As the AUC and TSS values of the 10 algorithms except SRE are all above 0.8 and 0.7, respectively, we excluded SRE in subsequence analysis ( Table 1). Results of paired Wilcoxon signed rank test indicate that full SDMs performed significantly better than climate SDMs in terms of both AUC (V = 466, p < .001) and TSS (V = 835, p < .001). Therefore, we hereafter reported the results regarding ensemble models and variable importance only for full SDMs.
The AUC and TSS value of the full ensemble model is 0.978 and 0.909, respectively, which is higher than that of the individual algorithms, indicating that our model was statistically robust and the predictive performance was near perfect (Allouche et al., 2006).
According to the ensemble model, environmental variables contributed differently to our models ( Table 2). The temperature annual range is the most influential variable, contributed 62.0% in the model, followed by annual mean temperature (10.4%), the percentage of grassland cover (9.4%), and precipitation of the driest month (8.1%). Response curves of the above four variables indicate that D. involucrate occurs mainly in areas with annual mean temperature ranging from approximately −1.6 to 19.2°C, temperature annual range ranging from approximately 21.3-31.5°C, precipitation of the driest month between about 0.1 and 59 mm, and the proportion of areas covered by grassland <0.59 (Figure 2). The remain eight variables contributed little to the distribution of D.

| Projected range shifts under different future scenarios
Overall, both climate SDMs and full SDMs predicted severe range contraction of the future distribution of D. involucrate, but full SDMs, by integrating climate and land-use change scenarios, predicted far less range contraction than climate SDMs, which only involve climate change scenarios (Figures 3 and 4). Specifically, the relative change

| The future effectiveness of current PAs
The current PAs only cover 6.8% (61,500 km 2 ) of the current suitable habitat of the D. involucrate. The relative change ratios of projected suitable habitat inside current PAs have a similar trend in predicted species range size both in the whole study area and in the current suitable range of D. involucrate ( Figure 5)

| DISCUSS ION
Accurately predicting the future distribution of species is crucial for understanding how species will response to global environmental change and evaluating the effectiveness of current PAs (Wiens et al., 2009) and has required that projections should cover multiple factors (e.g., climate and land use change) that are important drivers of species distributional changes (Brook et al., 2008;Travis, 2003).

ACK N OWLED G M ENTS
This study was supported by the National Natural Science

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Datasets used in this study are obtain from Long et al. (2021a), which is available online from the Dryad Digital Repository: https://doi.